%% Code for Section 6 of Acemoglu and Loebbing (2023): Automation and Polarization
% This code replicates the calibration of the model as described in Section
% 6.1
% It produces the matlab data files "calibration", "calibrationgrid",
% "calibrationlocalsearch", and "calibrationmw" in the folder
% Output/Calibration Results
% see online Appendix B.5 for an in-depth explanation of the numerical solution
% procedure of the model, the calibration, and the data sources

% specify path to base folder (Replication AL2025):
basefolder = fullfile('C:','Users','jonas.loebbing','Desktop','Replication AL2025');

% get baseline warning state
basewarnstate = warning;

%% create parameter structure
p = parameters();

%% set external parameters
% elasticity of substitution between tasks, from Humlum (2021)
p.lambda = 0.5;

%% set calibration targets
% note: all wage levels are measured in 2008 USD
% wage data from Output/Wage Data/wage.csv
% capital income share from Data/Capital Shares and User Costs/Acemoglu
% Loebbing Capital Shares and User Cost.xlsx
% minimum wage data from Data/Minimum Wage/Acemoglu Loebbing Minimum
% Wages.xlsx

% load wage data
wagedata = readtable(fullfile(basefolder,'Output','Wage Data','wage.csv'));

% extract variables
sdata = table2array(wagedata(:,1))./100;
logwdata = table2array(wagedata(:,2));
deltalogwdata = table2array(wagedata(:,3));

% find indices of 10, 30, 50, 90 percentiles
ind10 = find(sdata>=0.1,1);
ind30 = find(sdata>=0.3,1);
ind50 = find(sdata>=0.5,1);
ind90 = find(sdata>=0.9,1);

% calibration targets
logw5010data = logwdata(ind50)-logwdata(ind10);
logw9050data = logwdata(ind90)-logwdata(ind50);
alphakdata = 0.16; % income share of non-residential equipment and software, excl. all other capital
deltalogw3010data = deltalogwdata(ind30)-deltalogwdata(ind10);
deltalogw5030data = deltalogwdata(ind50)-deltalogwdata(ind30);
deltalogw9050data = deltalogwdata(ind90)-deltalogwdata(ind50);
deltaalphakdata = 0.01; % change in income share of non-residential equipment and software, excl. other capital, 1980-2016
logw50data = logwdata(ind50);

% collect targets
% for calibration without minimum wage
targetdata = [logw5010data; logw9050data; alphakdata; deltalogw3010data; deltalogw5030data; deltalogw9050data; deltaalphakdata];
% for calibration with minimum wage, now including 1980 median wage level
targetdatamw = [logw5010data; logw9050data; alphakdata; deltalogw3010data; deltalogw5030data; deltalogw9050data; deltaalphakdata; logw50data];

% choose weights for deviations from targets
% without minimum wage
weights = ones(numel(targetdata),1);
% with minimum wage, one more target so need one more dimension
weightsmw = ones(numel(targetdatamw),1);

% minimum wage levels 1980 and 2016, in 2008 dollars
mw = 7.1; 
mw16 = 7.46; 

%% grid search, no minimum wage
% construct parameter grid, bounds based on initial exploration of
% parameter space using knowledge of the economic model
mrange = [0.5 2.5 5]; % range for m
nrange = [0 -2.35 -4.7]; % range for n
arange = [1 7 14]; % range for a
brange = [-1 -300 -600]; % range for b
akrange = [1 110 220]; % range for ak
sqrange = [0.005 0.1]; % range for sq: q = min(qmax,qzero + sq*(qmax-qzero)), qmax s.t. q_2016<q_inf

[mgrid,ngrid,agrid,bgrid,akgrid,sqgrid] = ndgrid(mrange,nrange,arange,brange,akrange,sqrange);

% reshape grid to obtain a matrix where each row is a unique parameter
% combination
param = [reshape(mgrid,[],1) reshape(ngrid,[],1) reshape(agrid,[],1) reshape(bgrid,[],1) reshape(akgrid,[],1) reshape(sqgrid,[],1)];

% initiate matrix for model targets and exitflags
targetmodel = zeros(numel(mgrid),12);

% compute model targets and exitflags on the grid
parfor i = 1:numel(mgrid)

    warning('off','all');

    disp([num2str(i),' started']);
    
    % compute targets from model
    t = targetfun(param(i,:),0,0,p);
    targetmodel(i,:) = [t.targets' t.efs'];
    
    disp([num2str(i),' finished']);

end

% deviations from targets
dev = targetmodel(:,1:size(targetdata,1)) - repmat(targetdata',size(targetmodel,1),1);

% weighted sum of absolute deviations
sad = sum(abs(dev.*repmat(weights',size(dev,1),1)),2);

% find minimum of sum of absolute deviations
[sadmingrid,indmingrid] = min(sad);

% obtain minimizing parameter vector
parammingrid = param(indmingrid,:);

% save model results, parameter grid, weights and parammin
save(fullfile(basefolder,'Output','Calibration Results','calibrationgrid.mat'),'targetmodel','param','sad','weights','parammingrid','sadmingrid','indmingrid');

%% local minimization, no minimum wage
% load results from grid search if not already in workspace:
if exist('targetmodel','var') == 0
    load(fullfile(basefolder,'Output','Calibration Results','calibrationgrid.mat'));
end

% represent parameters as combination of two (extreme) vectors, mapping
% them to a common scale
% extreme vectors:
paramlow = min(0,2*parammingrid);
paramup = max(0,2*parammingrid);

% search for local minimum of sum of absolute deviations
fminoptions = optimset('Display','iter');
initialp = [0.5 0.5 0.5 0.5 0.5 0.5]; % using parammingrid as starting value
warning('off','all');
tic
[paramtransminloc,sadminloc,efloc] = fminsearch(@(paramtrans) sadfun(paramtrans,paramlow,paramup,weights,targetdata,0,0,p),initialp,fminoptions);
toc
% map parameters back to their original scale
paramminloc = (1-paramtransminloc).*paramlow + paramtransminloc.*paramup;

% obtain model targets at local minimum
targetmodelminloc = targetfun(paramminloc,0,0,p);
warning(basewarnstate);

% save results of local search
save(fullfile(basefolder,'Output','Calibration Results','calibrationlocalsearch.mat'),'paramminloc','sadminloc','efloc','targetmodelminloc','weights','parammingrid');

%% scale the model to match wage levels
% without minimum wage, increasing A and reducing q by the same number
% leaves all results unchanged except for wage levels and output

% load results of local search if not already in workspace
if exist('targetmodelminloc','var') == 0
    load(fullfile(basefolder,'Output','Calibration Results','calibrationlocalsearch.mat'));
end

% extract log median wage in 1980
logw50 = targetmodelminloc.targets(8);

% adjust A
A = p.A + logw50data - logw50;

% include new A in parameter vector
parammin = [paramminloc A];
% note that q is represented by sq (see above) in paramminloc
% when adjusting A and holding sq fixed, q is adjusted in parallel to A

% check median wage level
warning('off','all');
targetmodelmin = targetfun(parammin,0,0,p);
logw50 = targetmodelmin.targets(8);
warning(basewarnstate);

% sum of absolute deviations remains unchanged because wage level is
% matched exactly
sadmin = sadminloc;

% recall exit flag from local search
ef = efloc;

% save results
save(fullfile(basefolder,'Output','Calibration Results','calibration.mat'),'parammin','sadmin','ef','targetmodelmin','weights','targetmodelminloc');

%% local minimization with minimum wage

% load results of calibration without minimum wage if not already in
% workspace
if exist('parammin','var') == 0
    load(fullfile(basefolder,'Output','Calibration Results','calibration.mat'));
end

% use results from calibration without minimum wage as starting value for
% calibration with minimum wage
paramstart = parammin;

% represent parameters as combination of two extreme vectors again, for
% common scale
paramlow = min(0,2*paramstart);
paramup = max(0,2*paramstart);

% search for local minimum of sum of absolute deviations
fminoptions = optimset('Display','iter');
initialp = [0.5 0.5 0.5 0.5 0.5 0.5 0.5];
warning('off','all');
tic
[paramtransmin,sadminmw,efmw] = fminsearch(@(paramtrans) sadfun(paramtrans,paramlow,paramup,weightsmw,targetdatamw,mw,mw16,p),initialp,fminoptions);
toc
paramminmw = (1-paramtransmin).*paramlow + paramtransmin.*paramup;

% obtain model targets at local minimum
targetmodelminmw = targetfun(paramminmw,mw,mw16,p);
warning(basewarnstate);

% save results of local minimization with minimum wage
save(fullfile(basefolder,'Output','Calibration Results','calibrationmw.mat'),'paramminmw','sadminmw','efmw','targetmodelminmw','weightsmw','parammin');